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We report Monte Carlo simulations of a lattice-polymer model that can account for both polymer 
crystallization and liquid-liquid demixing in solutions of semiflexible homopolymers. In our model, 
neighboring polymer segments can have isotropic interactions that affect demixing, and anisotropic 
interactions that are responsible for freezing. However, our simulations show that the isotropic 
interactions also have a noticeable effect on the freezing curve, as do the anisotropic interactions on 
demixing. As the relative strength of the isotropic interactions is reduced, the liquid-liquid demixing 
transition disappears below the freezing curve. A simple, extended Flory-Huggins theory accounts 
quite well for the phase behavior observed in the simulations. 

I. INTRODUCTION 



Lattice models of polymer solutions are widely 
used because of their simplicity and computational 
convenience. 1-8 When modeling a polymer solution, the 
polymer chain occupies consecutive sites on the lattice, 
each site corresponding to the size of one chain unit, while 
the remaining sites correspond to solvent. 

The use of lattice models for polymer solutions dates 
back to the work of Meyer 1 . Flory 2 and Huggins 3 showed 
how, using a mean-field approximation, the lattice model 
yielded a powerful tool to predict the solution properties 
of flexible 9-11 and semi-flexible 7 polymers. Various re- 
finements to the Flory-Huggins (F-H) model have been 
proposed by a number of authors (see, e.g. Refs. 4 to 6). 
FH-style models can account for liquid-liquid (L-L) phase 
separations with an upper critical solution temperature 
(UCST) driven by the site-to-site mixing pair interac- 
tions in polymer solutions - however, they are ill suited 
to describe polymer crystallization, i.e. liquid-solid (L-S) 
phase transitions. This limitation is not due to any in- 
trinsic drawback of polymer lattice models as such, but to 
specific choice for the polymer interactions in the original 
F-H theory. In fact, the factors that lead to polymer crys- 
tallization, i.e. interactions that favor compact packing 
and stiffness of the polymer chains, can be accounted for 
in a lattice model, by introducing anisotropic interactions 
between adjacent polymer bonds. 8 Clearly, in real poly- 
mer solutions, both crystallization and phase separation 
can occur on cooling. While lattice models for polymer 
solutions can account for both types of phase transitions, 
most theoretical and simulation studies have focused on 
one transition or the other, and less attention has been 
paid to their interplay. Such interplay may change the 
pathway of a phase transition 12 ' 13 and hence determine 
the complex structure-property relationships of mixtures 
containing crystallizablc polymers, which has been the 
subject of much experimental research dating back to 
Richards. 14 
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FIG. 1. Schematic phase diagram of a binary mixture with 
a conventional monotectic triple point. 

When the L-S phase-transition curve intersects the L- 
L coexistence curve, both curves are terminated at the 
resulting triple point. Below the triple point, the fluid 
phase may phase-separate into a dilute solution and a 
dense crystalline phase, as depicted in Fig.l. This com- 
bination of L-L demixing and crystallization is often re- 
ferred to as "monotectic" behavior and has been observed 
in many experiments. 10 ' 15,16 The morphology of polymer 
crystallites appears to be sensitive to the result of ther- 
modynamic competition on cooling. 17 Special attention 
has been focused on the monotectic triple point. The 
kinetic competition between L-L demixing and crystal- 
lization on cooling in the vicinity of this triple point is 
an important issue for sol-gel transition and membrane 
preparation. 18-20 On cooling through the triple point, 
L-L phase separation is expected to occur before crys- 
tallization, though both phase transitions have the same 
equilibrium temperatures. 21 As a consequence, the den- 
sity modulation produced during the early stage of L-L 
demixing may be frozen by subsequent crystallization. 22 
Such frozen-in density modulations can be a practical 
way to control the metastable morphology of polymer 
gels and membranes through thermally induced pro- 
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cesses. Therefore, the ability to predict phase diagrams 
of the type shown in Fig. 1 could be of considerable 
practical importance. 

In this paper, we study the interplay of polymer crys- 
tallization and L-L demixing using both mean-field the- 
ories and Monte Carlo simulations of simple lattice mod- 
els. In particular, we pay attention to the shift of the 
crystallization and L-L demixing curves in the phase di- 
agrams due to this interplay. 

The remainder of this paper is organized as follows: 
after an introductory description of the simulation tech- 
niques, we compare the simulation results with the rel- 
evant theoretical predictions for the L-L phase separa- 
tion curve without prior disorder-order phase transition 
on cooling. Next, we discuss the simulations and mean- 
field calculation of the L-S curves and its thermodynamic 
competition with L-L demixing. 



II. SIMULATION TECHNIQUES 

In our Monte Carlo simulations, wc used a singlc- 
site-jumping micro-relaxation model with local sliding 
diffusion 23 to model the time evolution of self- and 
mutually-avoiding polymers in a cubic lattice with pe- 
riodic boundary conditions. In this model, monomer dis- 
placements are allowed along both the cubic axes and the 
(body and face) diagonals, so the coordination number 
of each site includes all the neighboring sites along the 
main axes and the diagonals, and isg = 6 + 8 + 12 = 26. 
The single-site jumping model with cither kink genera- 
tion or end-to-end sliding-diffusion was first proposed by 
Larson et al. 24 The kink-generation algorithm was subse- 
quently developed into the bond-fluctuation model. 25 ' 26 
A hybrid model combining kink generation and sliding 
diffusion into one mode of chain motion, was suggested 
by Lu and Yang. 27 The present hybrid model considers 
sliding-diffusion moves that are terminated by smooth- 
ing out the nearest kink conformation along the chain, 23 
in accord with de Gennes's picture of defect diffusion 
along the chain. 28 It has been verified that this model 
correctly reproduces both static and dynamic scalings of 
short polymers in the melt. 29 

In our simulations, we consider systems containing a 
number of 32-unit polymer chains. The polymers reside 
in a cubic box with 32 3 lattice sites. The polymer con- 
centration was varied by changing the number of poly- 
mers in the simulation box. Monte Carlo sampling was 
performed using the Metropolis method. Three ener- 
getic parameters were used to model the intra- and inter- 
molecular interactions of the polymers. The first param- 
eter E c measures the energy penalty associated with hav- 
ing two non-collinear consecutive bonds (a "kink" ) along 
the chain; it is a measure of the rigidity of chains. The 
second parameter E p measures the energy difference be- 
tween a pair of parallel and non-parallel polymer bonds 
in adjacent, non-bonded positions. A positive value of E p 



favors the compact packing of parallel chain molecules in 
a crystal. Finally, the parameter B describes the energy 
penalty for creating a monomer- solvent contact. The to- 
tal change in potential energy associated with a Monte 
Carlo trial move is 

AE _ E c Ac + E p Ap + BAm 

k^f ~ k^f 

^Ac + Apf? ID 

where Ac denotes the net change in the number of kinks, 
Ap is the change in the number of non-parallel adjacent 
bonds, and Am measures the change in the number of 
monomer-solvent contacts, ks is the Boltzmann constant 
and T is the temperature. As shown in Eq. 1, three di- 
mensionless parameters control the acceptance probabil- 
ity of Monte Carlo trial moves: B/E c is the term that 
dominates the L-L demixing temperature but no effect 
at all on the freezing of the pure polymer system. In 
contrast, E p /E c completely determines the freezing tem- 
perature of the pure polymer system, but it has only a 
slight effect on the demixing temperature. In fact, from 
Eq. 6 below, it follows that, the critical demixing tem- 
perature is approximately a factor of q higher in the case 
E p = and B / ' E c =/= than in the case where the values 
of B and E p are interchanged. In what follows, E c j (fc^T) 
is used as a measure of the (inverse) temperature of the 
system. If E c is much larger than B and E p , the poly- 
mer chains behave as almost rigid rods. In contrast, if 
E c = 0, the polymers are fully flexible. In what follows, 
we chose E v j E c — 1 as a value typical for semiflexible 
chains. The choice of the value of B/E c (and thereby 
the L-L demixing region) is discussed in the following 
sections. In our simulations, we lowered the temperature 
by increasing the value of E c / (ksT) from zero in steps of 
0.002. At each step, the total number of trial moves was 
500 MC cycles, where one Monte Carlo cycle (MC cycle) 
is defined as one trial move per monomer. The first 400 
MC cycles at each temperature were discarded for equi- 
libration, after which samples were taken once per MC 
cycle and averaged. This process corresponds to a slow 
cooling of the sample system. 

The most direct way to establish the equilibrium phase 
diagram of this model system would be to compute the 
free energy of all phases. Here, we follow a different route: 
we attempt to locate the equilibrium phase-transition 
temperatures during the dynamic cooling process. How- 
ever, rapid cooling may lead to a significant supercooling 
mainly due to the presence of a free-energy barrier for 
homogeneous nucleation. This is particularly true in di- 
lute solutions and small systems. In order to identify the 
correct equilibrium coexistence curves in a dynamic cool- 
ing scheme, supercooling should be eliminated as much 
as possible. To this end, we introduced one solid layer of 
terraced substrate formed by extended chains, as shown 
in Fig. 2A. These terraces can induce heterogeneous nu- 
cleation with a very small free-energy barrier. On such 
a large, terraced substrate, layer-by-layer crystal growth 
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can take place directly, thereby obviating the need for ho- 
mogeneous nucleation. In order to increase the accuracy 
of the method near the onset of the phase transition, we 
monitored the properties of the system during successive 
blocks of 500 MC cycles. If, during such a block, we found 
evidence for the onset of a phase transition, we kept the 
temperature constant for a number of subsequent blocks, 
until no further drift in the system properties was ob- 
served. 

On cooling, the degree of order in the sample system 
can be traced by the Flory "disorder" parameter, defined 
as the mean fraction of non-collinear connections of two 



consecutive bonds along the chains. On the cubic lattice, 
where 24 out of 25 directions for the connection to the 
next bond are non- collinear, the high-temperature limit 
of the disorder parameter is 0.96. The degree of demix- 
ing of the system can be monitored by tracing the value 
of a "mixing" parameter, defined as the mean fraction of 
the sites around a chain unit, that are occupied by sol- 
vent. Our estimates of the onsets of phase transitions are 
based on the averaged results of five independent cooling 
processes characterized by the same energy parameters, 
but different seeds for the random-number generation. 
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FIG. 2. Effect of a terraced substrate on the onset of crystallization upon cooling. The figures shown above were obtained 
for a solution of model polymers with length r = 32, at a volume fraction <j> = 0.0625, with E p /E c = 1 and B/E c — 0. (A) 
Snapshot of an athermal sample system containing one layer of terraced substrate formed by extended chains, that are not 
included in the polymer volume fraction. Viewing along the extended chains. (B) Disorder-parameter cooling curves for the 
sample systems with a terraced substrate on cooling (solid line) and under the absence of a seed on cooling (dashed line). The 
arrow indicates the onset of phase transition. (C) Step-size dependence of the onset of crystallization on cooling. (D) Finite-size 
scaling of the onset of crystallization on cooling for the sample systems with denoted concentrations. All error bars are smaller 
than the symbols. The segments are drawn as a guide to the eye. 
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As can be seen from Fig. 2B, the presence of a ter- 
raced substrate significantly decreases the kinetic delay 
on cooling for polymer crystallization from a dilute solu- 
tion. The onset of crystallization induced by the terraced 
substrate becomes insensitive to the number of steps on 
the substrate when this number is larger than 8, see Fig. 
2C. One might expect that more steps on the substrate 
would cause the substrate to adsorb more chains. The 
fact that the phase-transition temperature becomes in- 
sensitive to the number of steps (here, and in what fol- 
lows, we use 32 steps), suggests that pretransitional ad- 
sorption has a negligible effect on the apparent phase- 
transition temperature. In contrast, if no "template" is 
present, the onset of crystallization from a dilute solu- 
tion, depends strongly on the system size. This effect is 
probably due to the volume dependence of the homoge- 
neous nucleation rate. It can be completely eliminated by 
the introduction of a terraced substrate, as demonstrated 
in Fig. 2D. 

In the following sections, we first consider the case that 
Ep/ E c is zero and hence no crystallization can take place, 
while B/E c is large enough to induce L-L demixing on 
cooling. Next, we switch on E p /E c . This allows us to 
study a phase diagram that exhibits both L-L demixing 
and freezing. 



III. RESULTS AND DISCUSSION 
A. Liquid-liquid demixing without crystallization 

If both B/E c and E p /E c are zero, the model only 
takes excluded-volume interactions and the temperature 
dependence of chain flexibility into account. Even in 
this case, the polymer solution may exhibit a disorder- 
order phase transition on cooling. 8 ' 30 This transition is 
not, strictly speaking, a freezing transition but rather an 
isotropic-nematic phase transition: it is induced by the 
increase in anisotropic excluded volume interactions be- 
tween polymer chains, due to the increase in chain rigid- 
ity on cooling. 7,31 This transition has recently been stud- 
ied extensively by Weber et al. 32 

If we increase the value of B/E c while keeping E p /E c 
equal to zero, we should reach a point above which L- 
L demixing occurs prior to the isotropic- nematic phase 
transition on cooling. 

We focused our attention on the L-L demixing with val- 
ues of B/E c beyond that critical value, and kept track of 
the "mixing" parameter on cooling. As the dense liquid 
phase wets the terraced substrate, the onset temperature 
of L-L demixing induced by such a substrate should be a 
good approximation to the equilibrium phase separation 
temperature. A tentative binodal curve can thus be ob- 
tained in simulations to compare with the predictions of 
mean-field theories. 

Figure 3 shows the binodal curves for the sample sys- 
tems with E p /E c = and B/E c = 0.25. The binodal 



curve can be estimated from the condition of equal chem- 
ical potential of the coexisting phases, using the Eqn. 2, 
the F-H expression for the mixing free-energy. 



k B T 



= (1 - 0)ln(l - , 



-ln(0) + 0(1 - 4>) 



(g-2)B 



where <j) is polymer volume fraction, r is the chain length, 
and q the lattice-coordination number. As can be seen 
from Fig. 3, the theoretical predictions show a small but 
constant deviation from the simulation results. 

Yan et al. 33 have shown that a second-order lattice- 
cluster theory may provide a better description of the 
critical point of the binodal curve obtained in computer 
simulations. To second order, the mixing free-energy 
change per lattice site is 6 



A/, 



, = (1 - 0)ln(l - <£) + ^H4>) 
kb-L r 



1 



-qe<t> 2 + C Q + C 1 e + C 2 e 2 



(3) 



where e = 2S/(fesT). Explicit expressions for Co, C\ 
and C2 in terms of </>, q and r are given in Ref. 6. 
When we compare the predictions of the second-order 
lattice-cluster theory with our simulations, (dashed curve 
in Fig. 3) we find that this theory does not lead to better 
agreement with the simulation data, except perhaps at 
high polymer concentrations. It should be noted that, 
for very long polymer chains, the lattice cluster theory 
may predict more than one critical point. 34 Hence, the 
predictions of this theory should be viewed with some 
caution. 35 
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FIG. 3. Liquid-liquid coexistence curves (Td) of the sam- 
ple system with E p /E c = and B/E c = 0.25. The solid 
line is calculated from the classical Flory-Huggins free-energy 
expression for polymer solutions, and the dashed line is cal- 
culated from the second-order expansion of the mixing free 
energy in lattice-cluster theory. The triangles are the onsets 
of liquid-liquid demixing induced by a terraced substrate on 
cooling. The error bars are smaller than the symbols. 
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B. Polymer crystallization and its interplay with 
liquid-liquid demixing 



Starting the calculation from Eq. 5, the F-H expres- 
sion for the mixing free-energy change becomes 



When we set B/E c = and E p /E c = 1, L-L demixing 
is preempted by freezing. In fact, an estimate based on 
mean- field theory (Eqn. 5 below) indicates that, for these 
parameter values, the freezing temperature of the pure 
polymer is a factor four higher than the critical demixing 
temperature. We assume that the onset of crystallization 
induced by the terraced substrate, yields a good approx- 
imation for the equilibrium melting temperature. It is 
this temperature that we subsequently compare with the 
corresponding prediction of mean-field theory. 

The mean-field expression for the partition function of 
the disordered polymer solution is given by 8 : 



A/, 



rt\ ri2 £ 



> z (r-2)n 2e (l- 



)n 2 (r-l)n 2 rn 2 



(4) 



where 

z c =l + (q- 2)exp( 

Z p = CXp[-- i 2^(l - 



*=«P(-*^) 

n± denotes the number of sites occupied by the sol- 
vent, ri2 the number of chains, each having r units, and 
n = ni + rn2. We note that, in the above expression, we 
have corrected an error in the expression for the partition 
function given in Ref. 8. The corresponding expression 
for the free-energy density (i.e. the Hclmholtz free energy 
per lattice site) is 

ffi = (l-0)ln(l-0) + ^ 
kb-i r 



2(r-l)w 2 
qn 



\ e p 1 
I k B T\ 



ln(qr/2) 



D 



+(*-2)^ + (l-l/r) 



- (1 - 2/r) lnz c + (1 - 1/r) 
q-2 E p 



2 k B T J 



k B T 



We assume that the pure polymer crystal is in its fully 
ordered ground state and that the partition function of 
this state is equal to one. In a pure polymer system, 
melting takes place at the point where the free energies 
of the crystal and the melt cross. For polymer solutions, 
the freezing curve can be computed by imposing that the 
chemical potential of the polymers in crystal and solu- 
tion arc equal, i. e. /i c — /j° = /i s — /z , where fjr is 
the chemical potential of polymers in the ground state. 
As the free energy of the crystal phase is assumed to be 
equal to zero, the chemical potential of the polymers in 
that phase is also equal to zero. The chemical potential 
of the polymers in solution is /i s = dF s jdn^- Thus by 
solving the equation d\nZ s jdni = by iteration, we can 
obtain the equilibrium melting temperature. 



= (1 - #n(l - 4>) + -m(</>) 
k B l r 



+cf>(l-4>)(q-2) 



k B T q r k B T ! 



The binodal L-L curves can be separately estimated with- 
out the consideration of L-S curves. 
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FIG. 4. Liquid-liquid demixing curves (denoted as Td) and 
liquid-solid transition curves (denoted as T m ) for the sam- 
ple system with variable energy parameter settings (denoted 
as T(E P /E C , B/E c )). (A) Theoretical curves calculated from 
Eq. 4 with q e ff — q — 2 (Flory-Huggins approach). Note 
that changing E v /E c from 1 to 0, leads to a 10% decrease in 
Td- In contrast, lowering B / E c by 0.15 reduced Td by more 
than 50%. An arrow indicates the position of possible triple 
point; (B) Onsets of phase transitions induced by a terraced 
substrate on cooling. The error bars are smaller than the 
symbols, and the segments are drawn to guide the eye. 

In Figs. 4 A and 4B, we compare the mean-field pre- 
dictions for the phase diagram with the simulation data. 
In view of the simplicity of the mean-field theory, the 
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agreement between theory (without adjustable parame- 
ters) and the simulation data, is gratifying. 

According to Eq. 6, we should expect that a positive 
value of E p /E c will increase the L-L demixing tempera- 
ture. This is precisely the behavior observed in Fig. 4, 
where the L-L demixing curve of the sample system with 
B/E c = 0.25 shifts up when the value of E p /E c changes 
from zero to one. By carefully choosing the parameters, 
such as B/E c = 0.1, we can "tune" the relative strength 
of the tendencies to crystallize and to demix, and observe 
the intersection of the L-S and L-L curves. 

Although a change in the value of B/E c cannot change 
the freezing temperature of pure polymers, it can change 
the L-S coexistence curve of polymer solutions. The rea- 
son is that a poor solvent favors phase separation (be it 
L-L or L-S). 

However, in the simulations, we observed that the L-S 
curves cross not only at <fi = 1 but also at a second point 
near <f) = 0.73. This crossing point is not related to the 
presence of the terraced substrate, as it has also been ob- 
served in the absence of such a template. 8 Possibly, this 
failure of the simple mean- field theory is due to the rather 
naive way in which it accounts for the effective coordina- 
tion of monomers. We point out that, in our estimate, 
we have assumed that the effective coordination number 
is equal to q — 2. But, in more sophisticated theoretical 
descriptions, q e ff (as in Eq. 3) is, itself, concentration 
dependent. 



mers, Ah u is the heat of fusion per chain-unit. The pre- 
dictions of Eq. 7 for the melting-point depression upon 
dilution have been verified by several experimental mea- 
surements at both high and low concentration ends. 31 ' 32 
The linear relationship predicted by Eq. 7 does hold 
for those simulations where L-L demixing does not oc- 
cur (see Fig. 5). According to Eq. 7, the values of 
—q e ffB/ Ah u and E c /Ah u can be obtained respectively 
from the slope and the intercept of the freezing "line" . 
We found that —q e ffB/Ah u depends nearly linearly on 
B/E c x E c /Ah u . In addition, Ah u /E c varies linearly 
with B/E c . Assuming that both relations are, in fact, 
linear, we find: q eff = 54.0 and Ah u = 41. OB + 13.0.E C 
respectively. The latter result implies a microscopic cou- 
pling between L-L demixing and polymer crystallization, 
consistent but not identical with the previous study. 8 

In this paper, we have addressed the equilibrium freez- 
ing and demixing curves of lattice polymers. In subse- 
quent work, we shall address the effect of the interplay be- 
tween demixing and freezing on the kinetics of the phase 
transformation. 
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Flory has proposed a semi-empirical relationship be- 
tween the melting point and the concentration of poly- 
mers in solutions, 30 as given by 



1 

_L -ff. 



1 



k B 



(7) 



where is the equilibrium melting point of bulk poly- 
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